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Abstract 

Pegasus  Runway,  located  13  km  south  of  McMurdo  Station,  Antarctica,  on 
the  McMurdo  Ice  Shelf  (MIS),  is  constructed  out  of  snow  and  ice.  It  is  sus¬ 
ceptible  to  weakening  and  damage  caused  by  melting  and  to  reduction  in 
the  strength  caused  by  warm  weather  and  sunlight.  This  report  describes 
the  development  of  the  Pegasus  Runway  temperature  model.  It  begins  by 
quantitatively  describing  the  physical  properties  of  the  Pegasus  Runway 
snow  and  ice  and  the  physical  properties  of  the  MIS  directly  beneath  the 
runway.  The  temperature  model  is  based  on  a  one-dimensional  heat  con¬ 
duction  model  that  includes  the  penetration  and  absorption  of  solar  radia¬ 
tion  beneath  the  surface.  The  report  describes  the  methods  for  estimating 
the  sensible  heat,  latent  heat,  shortwave  radiation,  and  long-wave  radia¬ 
tion  surface  heat  fluxes  that  drive  the  model  and  presents  estimates  of  the 
constant-temperature  lower-boundary  condition  for  the  model.  A  novel 
approach  for  estimating  the  initial  vertical  temperature  profile  is  used.  We 
simulate  the  Pegasus  Runway  temperatures  for  three  austral  summer  sea¬ 
sons  (2011-12,  2012-13,  and  2013-14).  The  model  simulation  shows  good 
results  when  compared  to  in  situ  observations  of  the  runway  tempera¬ 
tures. 
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1  Introduction 

The  Pegasus  Runway,  intended  for  heavy  wheeled  aircraft,  was  construct¬ 
ed  of  snow  and  ice  in  the  austral  summer  of  1992-93.  It  is  located  13  km 
south  of  McMurdo  Station,  Antarctica,  on  the  McMurdo  Ice  Shelf  (MIS),  a 
lobe  on  the  extreme  west  of  the  Ross  Ice  Shelf.  The  area  was  selected  after 
extensive  study  (for  an  overview,  see  Blaisdell  et  al.  [1998]  and  also  Mellor 
[1988]  and  Mellor  and  Swithinbank  [1989])  and  has  unique  snow  cover 
and  glaciological  characteristics  (Stuart  and  Bull  1963;  Paige  1968; 
Swithinbank  1970)  that  are  favorable  for  a  permanent  airfield  (Klokov  and 
Diemand  1995).  The  area  where  the  runway  was  constructed  has  a  “thin, 
but  permanent  and  complete,  snow  cover”  (Blaisdell  and  Lang  1995).  The 
3000  by  90  m  runway  was  constructed  in  several  steps  to  remove  the 
rough  surface  ice,  to  fill  low  areas  and  subsurface  melt  holes,  and  to  grade 
the  runway  flat  and  smooth  (Blaisdell  and  Lang  1995).  During  subsequent 
seasons,  the  winter  snowfall  was  used  to  cover  the  ice  runway  and  to  pro¬ 
vide  a  high  albedo  protective  layer,  or  “snowcap,”  for  the  glacial  ice.  The 
snowcap  is  a  layer  of  snow  compacted  and  allowed  to  sinter  so  that  it  is 
strong  enough  to  support  wheeled  aircraft. 

The  Pegasus  Runway  has  been  used  each  season  since  1993.  It  has  served 
as  one  of  three  airfields  for  the  U.S.  Antarctic  Program  air  support  system 
at  McMurdo  Station,  which  also  includes  the  Sea  Ice  Runway,  located  on 
the  sea  ice  in  McMurdo  Sound,  and  Williams  Field,  located  on  the  snow- 
fields  east  of  Pegasus  on  the  Ross  Ice  Shelf.  In  most  years,  air  operations 
at  McMurdo  Station  have  used  the  Sea  Ice  Runway  early  in  the  season  with 
a  switch  over  to  Pegasus  and  Williams  Field  in  late  November  to  early  De¬ 
cember  when  the  sea  ice  is  no  longer  competent  enough  to  support  large 
cargo  planes  such  as  the  C-17.  Typically,  Pegasus  Runway  has  handled  in¬ 
tercontinental  flights  of  heavier  wheeled  aircraft  originating  from  Christ¬ 
church,  New  Zealand,  while  ski  equipped  aircraft  (e.g.,  LC-13OS)  have  op¬ 
erated  out  of  Williams  Field.  In  2009-10,  a  skiway  was  constructed  at 
Pegasus  (Figure  1);  and  the  Sea  Ice  Runway  and  Pegasus  were  operated 
sequentially  during  the  summer  with  both  wheeled  and  ski-equipped  air¬ 
craft  using  each  airfield.  During  this  time,  Williams  Field  was  minimally 
prepared  to  serve  as  an  emergency  divert  landing  site  for  ski-equipped  air¬ 
craft.  (See  Haehnel  et  al.  [2013]  for  a  more  thorough  review  of  air  opera¬ 
tions  at  McMurdo  Station).  A  critical  part  of  efficient  air  operations  is  the 
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ability  to  use  large  capacity  wheeled  aircraft  at  McMurdo,  such  as  the  C-17, 
A-319,  and  B-757,  to  bring  cargo  and  passengers  from  New  Zealand  and 
Australia  while  dedicating  the  slower  and  smaller  ski-equipped  LC-13OS 
for  transport  of  passengers  and  cargo  from  McMurdo  to  South  Pole  and  to 
camps  on  the  Antarctic  continent.  The  Sea  Ice  Runway  handles  these 
wheeled  flights  early  in  the  season  when  the  temperatures  are  cold.  Dur¬ 
ing  the  warm  part  of  the  season  (December  and  January),  these  wheeled 
flights  are  serviced  at  the  glacial  ice  runway  at  Pegasus. 


Figure  1.  Pegasus  Runway. 


CONTROL  TOWER 

ROAD  TO  McMURDO 

\ 

('  ") 

'  ’A 

r  \\\ 

V- 

PEGASUS  SKIWAY 

220FT  X 10000FT 

- 4- - - * 

*•  7 1 

PEGASUS  COMPACTED  SNOW  - ' 

v\ 

RUNWAY 

220FT  X  10000FT 

\\ 

2012  -  2013 

PEGASUS  SINGLE  AIRFILED 

:• 

SITE  PLAN 

r 

\S  \ 

N 

Given  that  the  Pegasus  Runway  is  constructed  out  of  snow  and  ice,  it  is 
susceptible  to  weakening  and  damage  caused  by  melting  and  to  reduction 
in  the  strength  of  the  snowcap  caused  by  warm  weather  and  solar  radia¬ 
tion.  The  critical  impact  of  solar  radiation  at  the  Pegasus  Runway  site  has 
long  been  recognized.  Paige  (1968)  describes  subsurface  melt  pools  that 
were  discovered  beneath  blue  glacier  ice  in  the  western  part  of  MIS  during 
a  reconnaissance  of  this  area  for  alternative  airfields  in  the  austral  summer 
of  1965-66.  He  stated  that  the  pools  were  caused  by  the  greenhouse  effect 
of  intense  solar  radiation,  low  albedo  of  the  blue  glacier  ice,  and  “heat  ab¬ 
sorption  of  dark  objects.”  Mellor  and  Swithinbank  (1989)  also  described 
melt  features  in  this  area. 
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Because  of  these  and  other  observations,  there  are  a  number  of  studies  of 
the  subsurface  temperature  regime  in  snow  and  ice  in  Antarctica  (see,  for 
example,  Schlatter  [1972],  Brandt  and  Warren  [1993],  and  Liston  et  al. 
[1999]).  These  studies’  main  focus  has  been  to  explain  subsurface  melting 
at  below-freezing  air  and  surface  temperatures.  In  summary,  these  studies 
found  that  subsurface  melting  occurred  in  only  blue  glacier  ice  that  was 
not  snow  covered.  In  snow-covered  areas,  the  maximum  of  the  vertical 
temperature  profile  could  occur  beneath  the  surface,  and  the  maximum 
was  limited  to  shallow  depths  with  little  or  no  subsurface  melting.  The  dif¬ 
ferences  between  the  conditions  found  in  snow  covered  areas  and  those  in 
areas  of  exposed  blue  ice  could  be  explained  “largely  by  radiative  and  heat- 
transfer  interactions  resulting  from  differences  in  albedo,  grain-size,  and 
density  between  the  two  mediums”  (Liston  et  al.  1999).  The  ability  of  a 
surface  snow  cover  to  prevent  subsurface  melting  was  exploited  at  the 
Pegasus  runway  where  it  was  found  “by  trial . . .  that  natural  melt  features 
could  be  prevented  at  the  Pegasus  site  by  completely  covering  exposed  ice 
surfaces  with  . . .  snow”  (Blaisdell  et  al.  1998).  In  short,  the  protective 
properties  of  snow  come  from  its  ability  to  reflect  most  of  the  downwelling 
solar  radiation  and  to  effectively  scatter  the  solar  radiation  that  does  pene¬ 
trate  the  surface. 

During  the  2012-13  and  2013-14  seasons,  Pegasus  Runway  experienced 
significant  reductions  in  the  strength  of  the  snow  and  ice  layers  that  com¬ 
posed  the  runway,  leading  to  extended  shutdowns.  Also,  during  short 
warming  trends  that  cause  temporary  reduction  in  the  strength  of  the  air¬ 
field,  operations  at  the  Pegasus  Runway  are  shifted  to  the  “night”  when  the 
sun  is  low  in  the  sky  and  the  ice  is  colder  and  stronger. 

These  disruptions  to  the  runway  operation  have  led  to  a  program  designed 
to  forecast  strength  reductions  of  the  snow  and  ice  layers  that  compose  the 
runway.  The  first  step  of  this  program  was  installing  additional  sensors  at 
the  Pegasus  Runway.  The  meteorological  station  at  Pegasus  Runway  was 
augmented  with  sensors  to  measure  the  downwelling  broadband 
shortwave  and  long-wave  radiation,  two  of  the  most  important  compo¬ 
nents  of  the  surface  heat  balance.  In  addition,  temperature  sensors  were 
placed  at  four  locations  in  the  Pegasus  Runway  to  measure  hourly  temper¬ 
atures  at  seven  levels  from  4  to  16  in.  of  depth.  These  measurements 
showed  the  runway  temperatures  were  highly  dynamic  with  large  daily 
and  seasonal  variations. 
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The  second  step  of  this  program  was  the  development  of  a  system  to  fore¬ 
cast  the  occurrence  of  strength  reduction  in  the  airfield.  The  two  major 
components  of  the  forecast  system  are  a  temperature  model  of  the  upper 
layers  of  the  Pegasus  Runway  and  a  model  of  the  runway  strength.  Sepa¬ 
rate  reports  will  describe  the  runway  strength  model  and  the  forecasting 
system.  The  present  work  describes  the  development  of  the  temperature 
model  of  the  upper  layers  of  the  Pegasus  Runway.  This  report  first  de¬ 
scribes  the  snow  and  ice  properties  of  Pegasus  Runway  and  the  MIS  di¬ 
rectly  below  the  runway.  We  next  describe  the  development  of  the  Pegasus 
Runway  temperature  model.  This  model  is  based  on  a  one-dimensional 
heat  conduction  model  and  includes  the  penetration  and  absorption  of  so¬ 
lar  radiation  beneath  the  surface.  We  then  describe  the  methods  for  esti¬ 
mating  the  heat  fluxes  that  drive  the  model.  These  are  sensible  heat,  latent 
heat,  shortwave  radiation,  and  long-wave  radiation  surface  heat  fluxes. 
Finally,  we  simulate  the  Pegasus  Runway  temperatures  for  three  austral 
summer  seasons  (2011-12,  2012-13,  and  2013-14). 
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2  Pegasus  Runway 

2.1  Overview 

Pegasus  Runway  is  constructed  on  the  MIS,  which  moves  in  a  generally 
northern  and  western  direction  (Swithinbank  1970)  and  periodically 
breaks  off  into  McMurdo  Sound.  Haehnel  et  al.  (2013)  report  the  move¬ 
ment  at  the  Pegasus  Airfield  to  be  about  30  m  per  year.  Drilling  in  De¬ 
cember  1991  found  the  ice  shelf  to  be  approximately  33  m  thick  at  the 
north  end  of  the  Pegasus  runway  with  4.8  m  of  freeboard  above  sea  level 
(Arcone  et  al.  1994).  The  bed  of  McMurdo  Sound  is  roughly  600  m  below 
sea  level  in  this  area. 

The  area  of  MIS  east  of  Pegasus  is  an  accumulation  zone  where  more  snow 
falls  in  a  typical  year  than  ablates.  The  annual  accumulation  gradually  de¬ 
creases  moving  westward  towards  Pegasus  Runway,  becoming  zero  at  the 
transition  from  the  accumulation  area  to  the  ablation  area  (Mellor  and 
Swithinbank  1989;  Swithinbank  1970).  This  transition  zone,  about  2  km 
wide  and  running  in  a  north-south  direction,  was  the  location  chosen  for 
the  Pegasus  Runway.  The  snow  cover  in  the  transition  zone  “varies  from 
local  patches  near  the  ablation  region  to  over  a  meter  where  it  phases  into 
the  accumulation  zone.  The  Pegasus  runway  was  sited  just  to  the  accumu¬ 
lation  side  of  where  the  snow  cover  ceases  to  be  patchy”  (Blaisdell  et  al. 
1998).  The  transition  zone  is  also  referred  to  as  the  zone  of  superimposed 
ice  due  to  the  extensive  melting  during  some  summers  and  subsequent  re¬ 
freezing  at  the  same  location  (Klokov  and  Diemand  1995);  the  Pegasus 
runway  is  located  on  superimposed  ice  (Blaisdell  et  al.  1998). 

2.2  Snow  and  ice  properties  of  Pegasus  Runway 

To  model  the  temperatures  in  the  Pegasus  runway,  it  is  necessary  to  have  a 
description  of  density,  thermal  conductivity,  specific  heat,  and  the  broad¬ 
band  shortwave  radiation  extinction  of  the  snow  and  ice  that  compose  and 
underlay  the  runway.  The  surface  albedo  must  be  known  as  well.  It  is  not 
necessary  that  the  numerical  model  extends  through  the  entire  vertical 
thickness  of  the  MIS  at  the  Pegasus  Runway,  only  that  the  model  extends 
to  a  depth  at  which  a  temperature  boundary  condition  can  be  set.  Typical¬ 
ly,  one  would  expect  the  temperature  to  be  relatively  constant  and  equal  to 
the  average  annual  air  temperature  at  a  depth  of  about  10-15  m  below  the 
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surface.  However,  the  MIS  thickness  is  finite  with  a  lower  boundary  tem¬ 
perature  determined  by  the  ice-water  interface  temperature,  a  tempera¬ 
ture  much  warmer  than  the  average  annual  air  temperature.  That  means 
that  it  is  necessary  to  consider  the  conditions  throughout  the  entire  MIS 
thickness  to  estimate  a  temperature  boundary  condition  beneath  Pegasus. 

Measurements  of  the  physical  properties  in  the  field  are  limited  to  snow 
and  ice  density  along  with  a  few  observations  of  broadband  shortwave  ra¬ 
diation  extinction.  We  estimate  the  other  properties  based  on  the  density 
and  the  ice  type.  The  direct  measurements  of  density  made  at  the  runway 
by  us  and  by  others  were  confined  to  the  top  5  m.  We  infer  the  density 
deeper  than  5  m  based  on  measurements  made  by  others  in  the  vicinity  of 
the  Pegasus  Runway.  The  vertical  profile  under  the  Pegasus  runway  is  ba¬ 
sically  superimposed  ice  underlain  by  firn.  Firn  is  the  intermediate  stage 
in  the  transition  from  snow  to  ice.  At  the  upper  surface  of  the  superim¬ 
posed  ice  is  the  thin,  artificially  maintained  snowcap.  The  underlying  firn 
has  two  distinct  layers:  an  upper  “normal”  (dry)  firn  layer  and,  beneath 
that,  a  layer  of  brine-saturated  firn  that  extends  downwards  to  the  bottom 
of  the  MIS.  In  all,  there  are  four  distinct  layers  of  snow  and  ice  at  the  run¬ 
way  between  the  surface  and  the  bottom  of  the  MIS:  snowcap,  superim¬ 
posed  ice,  firn,  and  brine-saturated  firn. 

As  stated  previously,  measurements  in  1991  found  the  ice  shelf  to  be  ap¬ 
proximately  33  m  thick  at  the  north  end  of  the  Pegasus  runway  with  4.8  m 
of  freeboard  above  sea  level.  This  thickness  is  consistent  with  the  thick¬ 
ness  map  by  McCrae  (1984)  and  the  recent  airborne  estimates  of  the  MIS 
thickness  by  Rack  et  al.  (2013).  In  general,  the  thickness  of  the  MIS  de¬ 
creases  with  proximity  to  its  northern  ocean  boundary.  Kovacs  et  al. 
(1982)  measured  thicknesses  of  roughly  10-15  m  within  a  kilometer  of  the 
edge  and  over  90  m  about  10  km  from  the  edge.  The  decrease  in  thickness 
towards  the  Ross  Sea  is  apparently  caused  by  ablation  due  to  the  intrusion 
of  warm  surface  water  beneath  the  MIS  during  the  austral  summer 
months  (Stern  et  al.  2013).  In  contrast,  Rack  et  al.  (2013)  detected  ice 
platelet  layers  below  the  ice  shelf  in  the  areas  of  greatest  thickness,  indi¬ 
cating  areas  of  basal  freezing  and  supercooled  water  emerging  from  below 
the  central  ice  shelf  cavity  into  McMurdo  Sound. 
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2.2.1  Snowcap 

2.2. 1.1  Density 

The  snowcap  on  the  runway  surface  is  approximately  27  cm  thick.  It  is 
regularly  compacted,  planed,  and  dragged  throughout  the  summer  months 
of  November  through  January  to  maintain  a  highly  reflective,  porous  sur¬ 
face.  We  made  measurements  of  the  density  of  the  snowcap  as  a  function 
of  depth  at  18  locations  along  the  length  and  width  of  the  runway  during 
18-27  December  2011.  These  measurements  extended  to  a  depth  of  40  cm 
and  included  both  the  snowcap  (top  10-12  cm)  and  underlying  ice.  Figure 
2  presents  this  data;  a  least-squares  fit  to  the  data  is 

A=(A.-A,)(l-«-")  +  A,  (D 


where 

ps  =  snowcap  density  (kg  m-3); 
po  =  535-05  (kg  m-3); 
pice  =  density  of  ice,  917  (kg  m-3); 
z  =  depth  (m); 
k  =  6.98  (m-1)- 


Figure  2.  Measured  density  of  the  snowcap. 
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Two  of  the  measurements  indicate  densities  greater  than  the  density  of  ice, 
917  kg  m~3.  It  is  likely  that  these  samples  contained  sediments  that  have  a 
density  greater  than  ice.  While  overall  sediment  content  of  the  snowcap  is 
not  known,  sediment  has  a  significant  presence  on  the  MIS;  the  far  west¬ 
ern  edge  of  the  MIS  is  heavily  covered  with  sediment  and  debris.  Rack  et 
al.  (2013)  estimate  sediments  accumulated  at  the  surface  and  present 
within  the  ice  shelf  can  add  up  to  a  layer  thickness  of  2  m  in  that  area. 
Towards  the  eastern  edge  of  the  MIS,  Dunbar  et  al.  (2009)  found  an  annu¬ 
al  accumulation  rate  of  0.80  g  m~2  year-1  of  aeolian  sediment  onto  the  ice 
shelf  surface.  They  attribute  this  to  strong  storms  that  carry  unconsolidat¬ 
ed  sediment  from  snow-free  areas  to  the  south,  such  as  Black  Island. 

2.2.1.2  Thermal  conductivity 

The  thermal  conductivity  of  heavily  compacted  snow  is  not  known  with 
any  certainty.  The  density  of  the  snowcap  ranges  from  550  to  856  kg  m-3. 
These  values  are  greater  than  the  samples  used  by  Sturm  et  al.  (1997)  to 
develop  their  thermal  conductivity  relationship  for  seasonal  snow.  (In 
fact,  they  limited  the  application  of  their  equation  to  snow  densities  less 
than  600  kg  m-3.)  To  account  for  the  possible  differences  between  the 
heat  transfer  properties  of  the  Pegasus  Runway  snow  cover  and  seasonal 
snow  and  to  ensure  that  thermal  conductivity  transitioned  smoothly  to 
that  of  ice  as  the  density  increased  above  600  kg  m-3,  we  propose  the  fol¬ 
lowing: 

for  ps  <  600  kg  m-3, 

ks  =0.138-1.01x10' ~3ps  +  3.233  x  10' ~6p2s,  (2) 

and  for  ps  >  600  kg  m-3, 

ks=m(Ps-600)  +  k6  00-  (3) 

where 

k  -  T 

m  =  ice  600  , 

Pice  “  600 

ks  =  the  thermal  conductivity  of  the  compacted  snow  (W  [m  °C]-1), 
/c6oo  =  the  thermal  conductivity  of  snow  with  density  600  kg  m-3 
calculated  using  equation  (2). 
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This  approach  provided  similar  results  to  the  average  of  the  lower  and  up¬ 
per  limits  of  thermal  conductivity  of  firn  provided  by  Paterson  (1994). 

2.2. 1.3  Heat  capacity 

We  set  the  specific  heat  of  the  snowcap,  per  unit  mass,  to  equal  that  of  ice, 
2114  J  (kg  °C)-1. 

2.2. 1.4  Surface  albedo 

The  albedo  of  the  snow  cover  at  the  Pegasus  Runway  site  was  measured 
(M.  Knuth*,  pers.  comm.)  in  15-minute  intervals  from  29  October  2010 
until  5  February  2011.  Two  broadband  pyranometers  (Eppley  Precision 
Spectral  Pyranometer)  were  mounted  on  a  mast  immediately  adjacent  to 
the  runway.  One  measured  downwelling  radiation  and  the  other 
upwelling  radiation  from  the  snow  surface.  The  ratio  of  the  upwelling  to 
the  downwelling  is  the  albedo.  Values  influenced  by  snow  or  frost  on  the 
upward  looking  pyranometers  were  discarded.  The  albedo  values  varied 
throughout  each  day  and  had  longer  trends  of  slight  decline  or  increase  of 
up  to  thirty  days  length.  The  daily  variation  correlated  well  with  the  solar 
angle  (Woolf  1968),  as  expected.  The  reasons  for  the  longer  trends  were 
not  investigated  but  probably  resulted  from  changing  snow  conditions. 
The  average  albedo  values  by  solar  angle  were  well  described  by  the  pa¬ 
rameterization  proposed  by  Dickinson  et  al.  (1986)  (as  described  in  Gard¬ 
ner  and  Sharp  [2010]): 


«,=«*<  6o°+max 


0,0.4 


1  -a 


02  >60° 


— 1 

yVl  +  2  but  J 


(4) 


where 

aez<6o°  =  the  average  albedo  when  the  solar  angle  is  less  than  6o°  (found 
to  be  0.76  in  the  present  case), 
b  =  a  fitting  constant  equal  to  0.1  in  the  present  case, 
ut  =  the  cosine  of  the  solar  angle. 

Figure  3  summarizes  the  measurements  and  equation  (4).  The 
pyranometers  are  rated  up  to  a  solar  angle  of  only  8o°,  and  we  ignored 


*  U.S.  Army  Cold  Regions  Research  and  Engineering  Laboratory,  Hanover,  NH. 
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values  measured  at  solar  angles  greater  than  this  (the  sun  is  low  on  the 
horizon). 

Following  the  examples  of  Henneman  and  Stefan  (1999)  and  Douville  et 
al.  (1995),  the  surface  albedo  is  also  reduced  at  a  time  constant  rate  when 
the  liquid  water  fraction  of  the  surface  layer  is  greater  than  zero.  An  albe¬ 
do  reduction  factor,  arf,  is  subtracted  from  at  to  arrive  at  the  effective  sur¬ 
face  albedo,  ate- 


a,e  =  «,  “  a,f 


(5) 


where 


arf  =  0.068^  (6) 

and  tmeit  is  the  total  time  in  days  from  the  start  of  the  time  period  with  the 
melt  fraction  of  the  surface  layer  greater  than  zero.  The  factor  0.068  was 
determined  empirically.  If  the  melt  fraction  returns  to  zero,  then  arf  is  set 
to  zero  for  the  entire  time  that  the  melt  fraction  is  zero.  However,  tmeit  is 
not  reset  to  zero. 

Figure  3.  Measured  and  modeled  snow  surface  albedo  with  solar  angle. 
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2.2. 1.5  Solar  radiation  extinction 

Blaisdell  et  al.  (1998)  measured  bulk  extinction  coefficients  for  snow  at  the 
Pegasus  site  as  30  m-1  for  snow  with  a  density  of  approximately  250  kg 
m~3,  20  m_1  for  snow  with  a  density  of  approximately  350  kg  m~3;  and 
10  m_1  for  snow  with  a  density  of  approximately  450  kg  m~3.  Bintanja  and 
van  den  Broeke  (1995)  recorded  two  measurements  for  snow  in  Antarctica: 

2.5  m_1  for  snow  with  a  density  of  850  kg  m~3  and  17.1  m_1  for  snow  with  a 
density  of  400  kg  m~3.  Based  on  these  measurements,  we  propose  the  fol¬ 
lowing  description  of  the  bulk  solar  radiation  extinction  coefficient,  ju 
(m-1),  as  a  function  of  the  snow  density: 

H  =  81.45e_0004ft  .  (7) 


Figure  4  shows  the  results. 


Figure  4.  Bulk  extinction  coefficient  for  Antarctic  snow. 
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2.2.2  Superimposed  ice 

2. 2. 2.1  Density 

The  superimposed  ice,  found  immediately  below  the  snowcap,  is  about 
4.7  m  thick.  In  January  1993,  Arcone  et  al.  (1994)  conducted  a  ground- 
penetrating  radar  (GPR)  survey  of  the  superimposed  ice,  profiling  the  en¬ 
tire  length  of  the  runway  along  the  western  side.  Solid  and  largely  bubble- 
free  ice  was  found  in  the  top  2  m  of  ice  in  agreement  with  the  ice  core  ob¬ 
servations  of  Blaisdell  et  al.  (1992).  Arcone  et  al.  (1994)  detected  a  num¬ 
ber  of  distinct  linear  horizons  up  to  30  m  long,  caused  by  density 
inhomogenities  resulting  from  of  air  bubbles  trapped  in  the  ice,  between  3 
and  7  m  depth.  As  part  of  the  survey,  a  5  m  ice  core  was  removed  from  the 
runway  and  the  ice  density  was  measured  at  0.10  m  intervals.  The  ice  den¬ 
sity  was  at  its  maximum  near  the  surface  with  a  density  very  close  to  pure 
freshwater  ice.  The  average  density  in  the  upper  4  m  was  856  kg  m~3  with 
a  standard  deviation  of  25.8  kg  m~3  (excluding  anomalous  layers).  Two 
thin  layers,  approximately  0.10  m  thick  at  1.9  m  and  3.7  m  depth,  were  ob¬ 
served  that  had  anomalously  low  densities  of  600  kg  m~3  or  less.  Arcone 
et  al.  (1994)  assumed  that  these  were  representative  of  density 
inhomogenities  that  caused  the  linear  horizons  seen  in  the  radar  profile. 
The  density  of  the  superimposed  ice  declined  abruptly  between  4.7  and  4.8 
m  depth  to  about  700  kg  m~3  at  5  m  depth.  This  represented  the  transition 
to  the  firn  layer  below. 

2. 2. 2. 2  Thermal  conductivity 

The  density  of  the  superimposed  ice  at  Pegasus  is  less  than  pure  ice  due  to 
the  presence  of  air  bubbles.  The  thermal  conductivity  of  ice  with  air  bub¬ 
bles  is  estimated  as  (Schwerdtfeger  [1963],  as  presented  by  Ashton  [1986]) 

,  _  i2kice  +  kc,  )  Pice  -  2  {Ke  ~  K  ){Pice  ~  Phi  )  ,  ,g. 

{2Ke  +  K  )  Pice  +{Ke  ~  K  )(Pice  ~  Phi  ) 


where 

km  =  the  conductivity  of  ice  with  air  bubbles  (W  [m  °C]“1); 
ka  =  the  thermal  conductivity  of  air,  0.025  (W  [m  °C]_1); 
pbi  =  the  density  of  ice  including  air  bubbles. 
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This  results  in  an  estimate  of  2.07  W  (m  °C)_1  for  the  thermal  conductivity 
of  the  superimposed  ice. 

2.2.2.3  Heat  capacity 

We  set  the  specific  heat  of  the  superimposed  ice,  per  unit  mass,  to  equal 
that  of  ice,  2114  J  (kg  °C)-1. 

2.2.2. 4  Solar  radiation  extinction 

Grenfell  and  Maykut  (1977)  measured  the  bulk  extinction  coefficient  of  the 
penetrating  fraction  in  freshwater  ice  to  be  1.5  m_1,  which  is  similar  to  the 
value  Liston  et  al.  (1999)  estimated,  and  is  the  value  used  here  for  the  su¬ 
perimposed  ice  layer. 

2.2.3  Firn 

2.2.3. 1  Density 

The  “normal”  (dry)  firn  layer  beneath  the  Pegasus  Runway  lies  between 
4.8  and  9  m  depth.  The  MIS  is  composed  predominately  of  firn.  Kovacset 
al.  (1982)  measured  the  vertical  distribution  of  firn  density  at  several  loca¬ 
tions  on  the  MIS  approximately  14  km  to  the  northeast  of  Pegasus.  Figure 
5  shows  the  results  of  their  measurements  as  a  function  of  depth.  Their 
measurements  show  the  firn  density  increasing  with  depth,  and  all  the 
measurements  at  all  locations  collapse  to  a  single  curve  when  plotted  with 
depth.  The  sharp  interface  where  they  encountered  the  brine-saturated 
firn  at  each  borehole  is  very  apparent  in  the  figure  where  the  density  ab¬ 
ruptly  increases  at  approximately  the  10,  20,  35,  and  50  m  depths.  Their 
density  measurements  in  the  dry  firn  are  well  described  as  a  function  of 
depth  by 

P firn  =  (  Pice  ~  Po  )  ( 1  “  )  +  Po  (9) 

where 

Pfim  =  firn  density  (kg  m-3); 
po  =  391-7  (kg  m-3); 
pice  =  density  of  ice,  917  (kg  m-3); 
z  =  depth  (m); 
k  =  0.0393  (m-1)- 
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Figure  5.  Observed  firn  density  of  MIS  with  depth  (Kovacs  et  al.  1982). 


However,  we  would  not  expect  the  density  beneath  the  Pegasus  Runway  to 
exactly  follow  the  density-depth  curve  shown  in  Figure  5.  The  overburden 
pressure  caused  by  the  superimposed  ice  on  the  underlying  firn  would  be 
greater  than  the  pressure  exerted  by  firn  of  the  same  thickness  and  vertical 
position  given  the  greater  density  of  the  superimposed  ice  compared  to 
firn.  Therefore  we  expect  the  overburden  pressure  to  be  a  better  indicator 
of  firn  density  beneath  the  runway  than  depth  beneath  the  surface  would 
be.  To  describe  the  firn  density  as  a  function  of  the  overburden  pressure, 
we  used  the  measurements  of  firn  density  by  Kovacs  et  al.  (1982)  and  es¬ 
timates  of  the  corresponding  overburden  pressure: 

P firn  =  (  Pice  ~  Po  )  (] 1  ~  ^  )  +  A)  (10) 


where 

po  =  444-52  (kg  m-3), 

P  =  overburden  pressure  (kPa), 
K  =  0.00575  (kPa-1) 
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Figure  6  shows  these  results;  note  that  the  data  points  influenced  by  brine 
have  not  been  included  with  this  graph. 


Figure  6.  Observed  firn  density  of  MIS  with  overburden  pressure  (Kovacs  et  al.  1982). 


The  overburden  pressure  at  the  bottom  of  the  superimposed  ice  is  roughly 

39.2  kPa  at  4.8  m  depth.  At  this  pressure,  one  would  expect  the  firn  densi¬ 
ty  to  be  540  kg  m-3  based  on  equation  (10).  According  to  equation  (9),  this 
density  would  be  found  at  a  depth  8.5  m  in  the  MIS.  We  estimated  the 
density  profile  of  the  firn  beneath  the  superimposed  ice  by  adding  3.7  m 
(the  difference  between  4.8  m  and  8.5  m)  to  the  depth  in  equation  (9)  for 
every  depth  below  the  superimposed  ice.  Based  on  the  GPR  measure¬ 
ments  of  Arcone  et  al.  (1994),  we  expect  the  transition  to  brine-saturated 
firn  to  occur  at  a  9  m  depth. 

2. 2. 3. 2  Thermal  properties 

We  assumed  that  the  functions  that  described  the  thermal  properties  of 
the  snowcap  also  apply  for  the  thermal  conductivity,  specific  heat,  and 
bulk  extinction  coefficient  of  the  normal  firn  layer. 
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2.2.4  Brine-saturated  firn 

The  layer  of  brine-saturated  firn  extends  from  a  depth  of  9  m  to  the  bot¬ 
tom  of  the  MIS  at  a  depth  of  33  m.  Brine-saturated  firn  was  first  detected 
in  the  MIS  in  1959  (Heine  1968).  The  brine  is  understood  to  originate 
from  seawater  entering  the  seaward  ice  front  of  the  MIS  after  breakout  of 
the  shelf  ice  into  McMurdo  Sound  (Heine  1968;  Risk  and  Hochstein  1967; 
Kovacs  and  Gow  1975;  Kovacs  et  al.  1982;  Cragin  et  al.  1983).  Apparently, 
permeable  ice  does  not  exist  at  the  bottom  of  the  ice  shelf;  so  brine  cannot 
enter  through  the  bottom  of  the  MIS.  As  a  result,  the  upper  surface  of  the 
brine  layer  is  consistently  at  a  lower  elevation  than  sea  level;  and  the  brine 
layer  surface  slopes  downwards  away  from  the  ice  front.  Using  GPR,  Ko¬ 
vacs  et  al.  (1982)  detected  step-like  waves  of  propagating  brine.  Each 
wave  corresponded  to  a  specific  breakout  event.  As  the  brine  moves 
through  the  firn,  it  becomes  more  concentrated  as  water  is  removed  by 
freezing,  finally  reaching  concentrations  of  five  and  seven  times  that  of  the 
original  seawater  (Cragin  et  al.  1983). 

Kovacs  et  al.  (1982)  mapped  the  limits  of  the  brine  penetration  in  the  MIS. 
This  map  indicates  that  the  location  of  Pegasus  Runway  is  likely  to  be  in 
the  region  of  brine-saturated  firn.  Arcone  et  al.  (1994)  also  confirmed  this 
in  their  January  1993  GPR  survey  of  the  Pegasus  runway  site.  The  domi¬ 
nant  feature  of  the  survey  was  a  strong  reflection  originating  between  8 
and  9  m  of  depth.  Arcone  et  al.  (1994)  attributed  this  layer  to  brine  that 
was  contiguous  with  the  brine  layer  identified  to  the  east  by  Kovacs  et  al. 
(1982).  We  look  to  the  measurements  of  Kovacs  et  al.  (1982)  for  infor¬ 
mation  on  the  densities  of  the  brine-saturated  firn  (Table  1).  First  of  all, 
the  measurements  show  an  abrupt  increase  of  density  at  the  interface  be¬ 
tween  the  dry  firn  and  the  brine-saturated  firn.  The  data  of  Heine  (as  re¬ 
ported  by  McCrae  1984)  also  shows  this  abrupt  transition.  Kovacs  et  al. 
(1982)  reported  measured  densities  for  three  cores  that  penetrated  the 
brine-saturated  firn  layer.  The  average  values  of  the  total  density  (brine 
and  ice)  ranged  from  836.6  to  877.8  kg  m~3  as  shown  in  Table  1.  Each  of 
the  cores  indicates  increasing  density  with  depth  within  the  brine- 
saturated  firn  layer.  All  the  cores  reach  a  maximum  of  about  880  kg  m~3 
1  m  below  the  top  of  the  brine  layer. 
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Table  1.  Measured  brine-saturated  firn  properties. 


Location 

Distance  from  MIS  edge(m) 

Brine  Depth  below  Surface  (m) 

Brine  Depth  below  Sea  Level  (m) 

MIS  Thickness  (m) 

Observed 

Calculated 

(Cox  and  Weeks  1983) 

Average  Melt  Salinity  (PPT) 

Average  Temperature  (C) 

Total  Density  (kg  rrr3) 

Ice  Fraction 

Brine  Fraction 

Brine  Salinity  (PPT) 

Air  Fraction 

Station  B 

800* 

8.85 

1 

20 

4.5 

-3.50 

877.8 

0.891 

0.058 

61.3 

0.051 

High  Side  of 
Brine  Step 

2840 

16.9 

4.4 

34 

14 

-14.72 

889.1 

0.896 

0.058 

174.5 

0.049 

Low  Side  of 
Brine  Step 

2900 

21.5 

8.8 

35 

25 

-12.75 

869.5 

0.811 

0.112 

162.0 

0.081 

Station  C 

3590 

19.4 

9.8 

40 

Station  D 

7370 

33.8 

13.2 

68 

8.3 

-15.0 

836.6 

0.891 

0.032 

177 

0.079 

Station  E 

9560 

50.4 

24.1 

84.6 

4.2 

-16.3 

871.5 

0.930 

0.015 

186 

0.056 

Station  F 

9700 

88.5 

*  estimated 


We  assume  an  abrupt  transition  to  brine-saturated  firn  at  a  depth  of  9  m 
beneath  Pegasus  runway  as  indicated  by  the  GPR  survey.  At  this  depth, 
the  density  jumps  abruptly  from  the  firn  density  to  that  of  the  brine- 
saturated  firn.  At  the  upper  surface  of  the  brine,  we  assume  a  total  density 
of  860  kg  m~3  with  a  linear  increase  over  the  next  meter  of  depth  to  a  den¬ 
sity  of  880  kg  m-3.  There  are  a  number  of  factors  to  take  into  account  to 
estimate  the  thermal  conductivity  and  specific  heat  of  the  brine-saturated 
firn.  First,  the  presence  of  brine,  ice,  and  air  all  impact  the  thermal  con¬ 
ductivity  and  heat  capacity  (Schwerdtfeger  1963;  Ono  1967;  Yen  et  al. 
1992).  Unfortunately,  we  have  no  measurements  of  the  actual  properties 
of  the  brine-saturated  firn  beneath  Pegasus  Runway.  We  expect  that  the 
density  of  the  original  dry  firn  prior  to  flooding  by  the  brine  would  follow 
the  pattern  described  by  equations  (9)  and  (10).  The  salinity  of  the  origi¬ 
nal  brine  that  flooded  this  location  is  not  known.  As  described  by  Cragin 
et  al.  (1983),  the  brine  can  become  more  concentrated  as  it  moves  through 
the  firn;  and  the  measurements  of  Kovacs  et  al.  (1982)  do  not  present  a 
consistent  picture  of  the  original  brine  salinity.  Freezing  of  the  brine 
would  form  additional  ice,  and  we  expect  that  the  fraction  of  ice  would  be 
greater  than  the  original  dry  firn  so  that  the  thermal  properties  would 
largely  reflect  those  of  ice.  However,  it  is  well  known  that  the  presence  of 
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brine  can  have  a  dramatic  impact,  especially  on  the  heat  capacity  at  warm 
temperatures  and  high  salinities.  We  approximate  the  thermal  properties 
based  on  the  properties  of  ice  with  bubbles  as  described  by  equation  (8). 
This  approximation  is  justified  given  the  9  m  depth  of  the  top  of  the  brine 
layer.  This  means  that  in  regions  where  the  temperature  is  changing  and 
heat  capacity  is  important,  the  ice  is  cold,  and  the  brine  has  less  impact  on 
the  thermal  properties.  Our  calculations,  discussed  in  Section  3.3  and  3.4, 
suggest  that  the  temperature  of  the  brine-saturated  firn  will  be  constant 
throughout  the  year  below  a  depth  of  15  m  beneath  the  surface  of  the  Peg¬ 
asus  Runway.  We  set  the  lower  boundary  of  our  simulation  domain  at 
15  m  where  we  can  assume  a  constant  temperature.  We  assumed  that  the 
thermal  conductivity  for  the  saturated  brine  layer  above  this  depth  was 
equivalent  to  that  of  ice  with  a  density  of  880  kg  m-3  and  calculated  the 
thermal  conductivity  as  2.13  W  (m  °C)  based  on  equation  (8). 

2.2.5  Summary  of  snow  and  ice  properties 

Figures  7  through  9  show  the  vertical  profiles  of  density,  thermal  conduc¬ 
tivity,  and  bulk  extinction  coefficients  used  in  the  Pegasus  Model.  The  ver¬ 
tical  density  profile  is  the  black  line  shown  in  Figure  7.  The  open  symbols 
are  the  data  points  used  to  develop  the  line  in  each  of  the  four  sections. 

The  dark  green  circles  are  the  measurements  of  Kovacs  et  al.  (1982)  of  the 
firn  density  near  the  Pegasus  Runway.  The  open  green  circles  are  the  data 
points  modified  to  account  for  the  additional  overburden  pressure  of  the 
superimposed  ice. 

The  vertical  thermal  conductivity  profile  is  the  black  line  shown  in  Figure 
8.  The  open  symbols  are  the  thermal  conductivities  calculated  for  the  den¬ 
sities  shown  in  Figure  7. 

The  vertical  bulk  extinction  coefficient  profile  is  the  black  line  shown  in 
Figure  9.  The  open  symbols  are  the  bulk  extinction  coefficients  calculated 
for  the  densities  shown  in  Figure  7.  We  expect  minimal  solar  radiation  to 
penetrate  below  3  m  of  depth  but  show  the  bulk  extinction  coefficient  pro¬ 
file  down  to  15  m  as  required  by  the  model. 
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Figure  7.  Observed  density  with  depth,  and  Pegasus  Model  density. 


Figure  8.  Estimated  thermal  conductivity  for  each  observed  density,  and  Pegasus  Model 

thermal  conductivity. 
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Figure  9.  Estimated  bulk  extinction  coefficient  for  each  observed  density,  and  Pegasus  Model 

bulk  extinction  coefficient. 
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3  Temperature  Estimation 


There  are  several  components  required  for  estimating  the  vertical  temper¬ 
ature  distribution  in  the  Pegasus  Runway.  The  first  component  is  the  heat 
transfer  model  of  the  runway  snowcap,  superimposed  ice,  firn,  and  brine- 
saturated  firn.  The  heat  transfer  model  is  based  on  a  one-dimensional 
heat-conduction  equation  that  includes  internal  heat  generation  through 
solar  radiation  absorption.  Particular  attention  must  be  placed  on  the  heat 
conduction  parameters  of  density,  thermal  conductivity,  and  radiation  ab¬ 
sorption,  which  vary  significantly  with  depth.  The  second  component  is 
estimating  the  model  upper  and  lower  boundary  conditions.  The  energy 
flux  absorbed  at  the  surface  forms  the  upper  boundary  condition  of  the 
heat  conduction  model.  The  lower  boundary  condition  is  a  fixed  tempera¬ 
ture.  A  third  component  is  the  initial  temperature  profile  that  must  be  set 
at  the  start  of  each  simulation  period.  We  discuss  each  of  these  compo¬ 
nents  below. 


3.1  Heat  conduction  model 

The  temporal  evolution  of  the  vertical  temperature  distribution  is  de¬ 
scribed  by  the  one-dimensional  heat  transfer  equation: 


P  c 

‘  "  ft  dz 


ks%_ 

dz 


dz 


(11) 


where 

Ti  =  the  temperature  of  the  fth  snow  or  ice  layer  (°C), 
z  =  the  depth  (m), 
t  =  time  (s), 

pi  =  the  density  of  the  ith  snow  or  ice  layer  (kg  m~3), 

Cp  =  the  specific  heat  (J  [kg  °C]“1), 
ki  =  the  thermal  conductivity  (W  [m  °C]_1), 

Qsw  =  the  solar  radiative  flux  (W  m-2). 

This  heat  conduction  equation  is  solved  using  a  Crank-Nicholson  finite- 
difference  scheme  following  the  approach  of  Flato  and  Brown  (1996).  This 
approach  allows  for  a  variable  solution  grid  that  divides  the  snow  and  ice 
composing  the  runway  and  the  MIS  beneath  the  runway  into  layers  of  ar- 
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bitrary  thickness.  Each  layer  is  assumed  to  have  a  uniform  thermal  con¬ 
ductivity  and  heat  capacity.  The  model  domain  extends  from  the  surface 
of  the  Pegasus  Runway  down  to  a  depth  of  15  m.  The  physical  properties 
of  the  runway  snow  and  ice  are  assumed  known  and  summarized  in  Fig¬ 
ures  7-9.  The  solution  requires  boundary  conditions  at  the  upper  surface 
and  at  the  bottom  surface  of  the  model  domain  and  a  known  initial  tem¬ 
perature  distribution  at  the  start  of  the  simulation  time  period.  The  upper 
boundary  condition  is  the  net  surface  heat  flux  absorbed  at  the  surface. 

The  lower  boundary  condition  is  a  constant,  known  temperature. 

3.1.1  Solar  radiative  flux 

It  is  well  known  that  the  albedo  and  absorption  of  solar  radiation  by  snow 
and  ice  varies  with  wavelength.  However,  given  the  operational  con¬ 
straints  at  the  Pegasus  runway,  only  broadband  solar  radiative 
downwelling  flux  data  is  currently  available  or  likely  to  be  available  in  the 
future.  The  shortwave  radiation  absorption  can  be  divided  into  two  pro¬ 
cesses:  the  absorption  at  the  surface  of  the  non-penetrating  fraction  and 
the  gradual  absorption  with  depth  of  the  penetrating  fraction  (Grenfell  and 
Maykut  1977;  Maykut  1982).  At  Pegasus  Runway,  the  surface  layer  is 
snow;  and  the  fraction  of  the  shortwave  radiation  penetrating  into  snow  is 
generally  small.  The  penetrating  fraction  of  the  shortwave  radiation,  Qsw  0, 
can  be  estimated  as 


Qswo=P(}-ate)Qswi 


(12) 


where 

P  =  the  penetrating  fraction, 

Qswi  =  the  downwelling  shortwave  radiation  (W  nm2), 
ate  =  the  effective  broadband  surface  albedo. 

The  albedo,  ate,  is  a  function  of  the  solar  angle  as  shown  in  equations  (4) 
and  (5)  and  is  reduced  by  the  presence  of  liquid  water  on  the  surface  layer. 
The  penetrating  fraction  is  set  at  0.17  following  the  example  of  Flato  and 
Brown  (1996)  for  thin  snow. 

The  absorption  of  the  penetrating  fraction  is  described  by  the  bulk  extinc¬ 
tion  coefficient,  which  varies  with  depth  as  shown  in  Figure  9.  The  broad- 
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band  absorption  of  the  penetrating  fraction  of  solar  radiation  absorption 
with  depth  is  described  by 


Q  =Q  ne~M:  (13) 

z-'SW  z  z-'SW  0  '  ' 

where  Qswz  is  the  flux  of  shortwave  solar  radiation  at  depth  z  in  the  snow- 
pack  (W  m-2)  and  //  is  the  bulk  extinction  coefficientCm-1). 

3.1.2  Impact  of  meltwater  on  physical  properties 

The  three  physical  properties  of  the  snow  and  ice,  pi,  Cp,  and  ki,  are  all 
modified,  following  the  example  of  Liston  et  al.  (1999),  when  meltwater  is 
present.  The  thermal  conductivity  is  modified  according  to 

kimel,=(l~fw)ki+fwK  (14) 


where 

fw  =  the  water  fraction  ( 1  >  fw>  0 ), 

ki  =  the  thermal  conductivity  of  layer  with  no  melt  (W  [m  °C]_1), 
kw  =  the  thermal  conductivity  of  water  (W  [m  °C]_1). 

A  similar  procedure  is  followed  for  pi  and  Cp. 

3.2  Net  surface  heat  flux 

3.2.1  Overview 

The  net  heat  flux  absorbed  at  the  surface,  F0,  serves  as  the  upper  boundary 
condition  for  the  heat  conduction  model 


td-L 

dz 


=  Fn 


surface 


(15) 


F0  is  composed  of  four  modes  of  heat  transfer:  the  net  long-wave  and 
shortwave  radiation;  the  latent  heat  flux,  Hv,  and  the  sensible  heat,  Hs. 


F.  =  Q,wl  -  ear;  +  (1  -  «.)  (1  -  fi)  +  HL  +  Hs 


(16) 
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where 

Qiwi  -  the  down  welling  long -wave  radiation  (W  m-2), 

Ts  =  the  surface  temperature  (°C), 
e  =  the  snow  surface  emissivity, 
cr  =  the  Stefan-Boltzmann  constant  (W  m~2  K_4), 
at  =  the  broadband  shortwave  surface  albedo, 

B  =  the  penetrating  fraction  of  the  shortwave  radiation, 

Qswi  -  the  down  welling  shortwave  radiation  (W  m-2). 

Both  Qiwi  and  Qswi  were  measured  at  the  Pegasus  Runway. 

3.2.2  Sensible  and  latent  heat  fluxes 

Modern  methods  for  estimating  surface  fluxes  of  sensible  heat,  Hs,  and  la¬ 
tent  heat,  Hl,  from  extensive,  horizontally  homogeneous  snow  and  ice  sur¬ 
faces  rely  on  a  bulk  flux  algorithm  (see  Andreas  [1996]  for  a  complete  re¬ 
view).  The  bulk  flux  algorithm  must  include  predictions  for  all  three 
turbulent  fluxes— momentum,  r,  and  sensible  and  latent  heat— because  the 
prognostic  equations  are  coupled.  The  basic  equations  are 


r  =  pCDrS;, 

(17) 

H  =  PaCpaC„rSr(Qs  -@r), 

(18) 

Hl  =  paLvCErS,.(Qs-Qr) 

(19) 

where 


Sr 

0r and  Qr 
0s  and  Qs 

Pa,  Cpa,  and  Lv 


wind  speed  at  reference  height  r  (m  s_1); 
the  potential  temperature  (°C)  and  specific  humidity 
(kg  kg-1)  of  the  air  at  r; 

the  temperature  (°C)  and  specific  humidity  (kg  kg-1)  of 
the  air  at  the  surface; 

the  air  density  (kg  m-3),  specific  heat  of  air  (J  [kg  °C]-1), 
and  latent  heat  of  vaporization  or  sublimation  (J  kg-1). 


The  temperature  of  the  air  at  the  surface  is  assumed  to  be  equal  to  the  sur¬ 
face  temperature  calculated  by  the  heat  conduction  model.  The  specific 
humidity  of  the  air  at  the  surface  is  determined  by  the  saturation  vapor 
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pressure  over  ice  at  the  calculated  surface  temperature  and  observed  air 
pressure  (Andreas  2005). 


The  key  to  using  equations  (17),  (18),  and  (19)  is  estimating  the  drag  coef¬ 
ficient  appropriate  at  reference  height  r  (CA-)  and  the  transfer  coefficients 
for  sensible  (Chv)  and  latent  (C&)  heat  at  r.  These  are  obtained  by  estimat¬ 
ing  the  roughness  lengths  for  wind  speed  (z0),  temperature  (zr),  and  hu¬ 
midity  ( zQ )  (e.g.,  Andreas  1998): 


Cor  = 


[ln(r/z„)  —  y/m  (r/Z)] 


(20) 


CHr  = 


[ln(r/z0)-^m(r/Z)][ln(r/zr)-^(r/Z)]  ’ 


(21) 


Cet  ~ 


[in  (r/z0)-^m(r/Z)][ln(r/ze)-^(r/Z,) 


(22) 


Here, 


k  =  0.4,  the  von  Karman  constant; 

L  =  the  Obukhov  length,  a  stratification  parameter  (m); 
y/m  and  y/h  =  known  functions  of  r/L. 

For  y/m  and  y/h,  we  use  Paulson’s  (1970)  functions  for  unstable  stratifica¬ 
tion  (i.e.,  for  L  <  o)  and  Holtslag  and  De  Bruin’s  (1988)  (see  also  Andreas 
[1998])  for  stable  stratification  (i.e.,  for  L  >  o). 

The  Obukhov  length  is  a  stratification  parameter;  in  the  algorithm,  we 
compute  it  as 


L  =  - 


O 

kg 


0.610 

/*  "i"  —  Cf  * 

1  +  0.610 


(23) 


V 
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Here, 

u*  =  friction  velocity  (m  s_1); 

f*  =  surface  layer  temperature  scale, - - —  (°C); 

p  C  u * 

r  a  pa  * 

q*  =  surface  layer  humidity  scale, - L —  (kg  kg'1); 

PaLvu* 

jj  =  gravity  (m  s~2); 

©  =  the  surface  layer  mean  temperature  (°C); 

Q  =  the  surface  layer  mean  humidity  (kg  kg-1)  • 

For  the  snow-  and  ice-covered  area  of  the  Pegasus  site,  we  use  a  constant 
value  for  z0,  l  mm.  This  is  a  typical  value  for  snow-covered  ice  (e.g.,  An¬ 
dreas  and  Claffey  1995).  For  zr  and  zQ  in  equations  (21)  and  (22),  we  use 

Andreas’s  (1987)  parameterization.  This,  likewise,  has  proven  useful  for 
estimating  heat  and  vapor  fluxes  over  surfaces  of  ice  or  snow  (e.g.,  Jordan 
et  al.  1999,  2001). 

The  Liu  et  al.  (1979)  and  Andreas  (1987)  models  are  somewhat  similar  in 
their  approaches  to  predicting  the  roughness  lengths  for  temperature  and 
moisture,  zT  and  zQ.  Both  of  these  define  a  roughness  Reynolds  number  as 

Rt  =  (24) 

v 

This  is  formed  like  the  usual  Reynolds  number— a  velocity  scale  times  a 
length  scale  divided  by  the  kinematic  viscosity  of  the  air,  v.  Here,  howev¬ 
er,  the  velocity  and  length  scales  are  related  to  the  friction  velocity  and  z0. 

The  Andreas  (1987)  model  predicts  the  scalar  roughness  length,  zs  (either 
zror  Zq),  from 


ln(zs/z0)  =  b0  +Z>jlni?*  +  A(lnR,)  ,  (25) 

where  Table  2  gives  the  polynomial  coefficients  (also  tabulated  in  Andreas 
[1987,  2002]). 
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Table  2.  Values  of  the  coefficients  to  use  in  equation  (25). 


Coefficient 

/?*<  0.135 

0.135  </?*<  2.5 

2.5  </?*<  1000 

Smooth 

Transition 

Rough 

Temperature  (zt/zo) 

bo 

1.250 

0.149 

0.317 

bi 

0 

-0.550 

-0.565 

b2 

0 

0 

-0.183 

Humidity  (zq/zo) 

bo 

1.610 

0.351 

0.396 

bi 

0 

-0.625 

-0.512 

b2 

0 

0 

-0.180 

Finally,  we  need  to  discuss  the  wind  speed  variable,  Sr,  in  equations  (17), 
(18),  and  (19).  When  the  stratification  is  unstable,  we  model  Sr  by  using 
the  COARE  gustiness  parameterization  (Fairall  et  al.  1996): 

S2=U2  +  /32w2.  (26) 


Here, 

Ur  =  the  wind  speed  measured  at  reference  height  r  (m  s_1), 
in*  =  the  convective  velocity  scale  (m  s_1)  (e.g.,  Deardorff  1970; 

Fairall  et  al.  1996), 

(3  =  1.25. 

Equation  (26)  acknowledges  that,  in  unstable  stratification,  convectively 
driven  atmospheric  motions  still  cause  turbulent  transport  even  when  the 
mean  winds  are  light. 

Similarly,  in  stable  stratification,  we  model  Sr  as 

5,.  =  Ur  +  W0 ,  (27) 

where  W0  is  a  so-called  windless  parameter  that  we  take  as  0.5  m  s_1  (Jor¬ 
dan  et  al.  1999).  Like  equation  (26),  equation  (27)  acknowledges  that, 
even  when  the  mean  wind  speed  approaches  zero,  in  stable  stratification 
the  atmosphere  still  has  ways  to  carry  out  vertical  exchange.  In  a  stably 
stratified  atmosphere,  gravity  waves  are  probably  the  main  agents  for  this 
“windless”  exchange. 
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Many  of  the  equations  in  our  algorithm  are  coupled  because  they  depend 
on  u*.  The  Obukhov  length  L  that  appears  in  (23)  is  the  main  coupling 
variable,  though,  because  it  contains  u*,  Hs,  and  Hl-  Thus,  we  must  solve 

this  system  of  equations  iteratively  by  first  assuming  neutral  stratification 
(i.e.,  r/L  =  o)  and  then  iteratively  calculating  the  actual  value  of  L. 


3.3  Lower  thermal  boundary  condition 


The  lower  thermal  boundary  condition  is  set  as  a  constant  temperature  at 
the  maximum  depth  of  the  solution  domain.  This  is  the  minimum  depth 
at  which  the  ice  temperature  is  constant  throughout  the  simulation  period. 
There  are  no  measurements  of  the  temperature  deep  beneath  the  runway 
to  use  as  a  guide  for  selecting  the  depth.  To  estimate  the  minimum  depth 
with  constant  temperature,  we  solved  equation  (11)  for  the  temperature 
profile  through  the  entire  thickness  of  the  MIS  beneath  the  Pegasus  Run¬ 
way.  This  simulation  used  a  simplified  surface  heat  flux  calculation  that 
assumed  that  the  surface  temperature  was  equal  to  the  air  temperature 
and  that  the  air  temperature  throughout  each  year  could  be  estimated  us¬ 
ing  cosine  functions.  First,  to  determine  the  cosine  functions,  the  average 
daily  temperature  of  each  day  of  the  year  at  the  McMurdo  station  was  cal¬ 
culated  based  on  the  period  of  record  (1973  to  the  present).  Then,  the  an¬ 
nual  cycle  of  temperatures  was  estimated  using  cosine  functions.  We 
found  that  only  two  periods  were  necessary  to  provide  a  very  reasonable  fit 
to  the  data  as  shown  in  Figure  10.  The  average  temperature  on  day  j,  Tj,  is 
found  as 


T.  =T+YjA> 


cos 


r  n2n  j 
v  365 


(28) 


where 


j  =  the  Julian  day  of  the  year, 

T  =  the  annual  average  temperature  at  McMurdo  Station 
(-17.04°C), 

An  =  the  amplitude  of  the  nth  cycle  (Ai  =  n.i9°C ,A2  =  4.29°C), 
£n  =  the  phase  of  the  nth  cycle  (ei  =  0.08,  e2  =  -0.27). 
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Figure  10.  Observed  annual  cycle  of  daily  average  temperatures  at  McMurdo  Station  with 

fitted  period  estimation. 


The  vertical  temperature  profile  was  estimated  using  one-day  time  steps. 
The  physical  characteristics  used  were  as  described  above,  and  the  condi¬ 
tions  at  15  m  were  extended  to  the  bottom  of  the  MIS  at  34  m  depth.  The 
thermal  boundary  condition  at  the  bottom  of  the  MIS  was  fixed  at  the 
ice/sea  water  equilibrium  temperature  of  -i.7°C.  The  simulation  started 
with  the  initial  vertical  profile  temperature  set  at  the  annual  average  air 
temperature  of  -17.04°C.  After  simulating  30  years,  during  which  the  sur¬ 
face  temperature  was  varied  according  to  equation  (28),  the  results 
showed  that  the  temperature  in  the  upper  portion  of  the  MIS  varied  daily 
in  response  to  the  variation  in  air  temperature;  and  the  temperature  on 
each  day  of  the  year  was  identical  to  the  year  prior.  The  temperatures  be¬ 
low  15  m  were  essentially  constant  with  time.  Figure  11  shows  the  results 
of  this  simulation,  displaying  the  temperature  profile  on  the  first  of  every 
other  month.  A  second  simulation  was  run  to  determine  a  steady  tem¬ 
perature  profile  through  the  entire  thickness  of  the  MIS  beneath  the  Pega¬ 
sus  Runway.  The  simulation  started  with  the  initial  vertical  profile  tem¬ 
perature  set  at  the  annual  average  air  temperature  of  -17.04°C.  The 
surface  temperature  was  fixed  at  the  annual  average  temperature  at 
McMurdo  Station  (-17.04°C),  and  the  MIS  bottom  temperature  was  fixed 
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at  a  constant  -i.y°C.  The  simulation  was  propagated  through  time  until 
steady  state  was  reached.  Figure  n  also  shows  the  results  of  this  second 
simulation.  The  overall  results  indicate  that  at  a  15  m  depth,  the  tempera¬ 
ture  remains  at  a  near  constant  temperature  of  -9-42°C  and  that  the  tem¬ 
peratures  below  this  depth  are  steady  and  equal  to  the  long-term  steady 
temperature  profile  based  on  the  annual  average  temperature  at  McMurdo 
Station  and  known  bottom  temperature. 


Figure  11.  Estimated  temperature  profile  throughout  the  MIS  beneath  the  Pegasus  Runway. 
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3.4  Initial  temperature  profile 

We  used  a  novel  approach  to  estimate  the  initial  temperature  profile 
through  the  Pegasus  runway  down  to  the  15  m  depth.  This  approach  is 
based  on  the  analytic  determination  of  the  propagation  of  periodic  oscilla¬ 
tions  of  the  surface  temperature  into  the  depth  first  given  mathematical 
form  by  William  Thomson  (described  in  Carslaw  and  Jaeger  [1980,  81]) 
combined  with  knowledge  of  the  long-term  steady  temperature  profile 
shown  in  Figure  11. 


The  approximate  temperature  at  depth  z  on  day  j  is  then 


ERDC/CRREL  TR-15-2 


31 


rjiZ  rj- r 


+  TS  +  Z  4 


-kzv« _ 

e  cos 


/?2;r  j 
365 


(29) 


where  is  the  steady  state  temperature  at  depth  z  (shown  in  Figure  li), 
Ts  is  the  mean  temperature  difference  between  McMurdo  Station  and  Peg¬ 
asus  Runway;  and 


I 


2  n  p£ 


k  = 


n\365  2k. 


i  J 


is  a  depth-averaged  value  of  the  thermal  properties  of  the  m  layers  from 
the  surface  to  the  depth  z. 

The  results  using  this  approach  usually  agreed  to  within  i°C-2°C  of  the 
long-term  simulations  described  in  the  previous  section. 
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4  Operational  Field  Data 

Meteorological  data  collected  at  the  Pegasus  Runway  was  used  as  input  to 
the  Pegasus  Runway  vertical  temperature  model.  The  meteorological  data 
included  wind  speed  (m  s_1),  air  temperature  (°C),  relative  humidity  (%), 
air  pressure  (mbar),  down  welling  solar  radiation  (W  m~2),  and 
downwelling  long-wave  radiation  (W  m~2).  The  U.S.  Navy  SPAWAR  (Space 
and  Naval  Warfare  Systems  Command)  Systems  Center  Atlantic,  Polar 
Programs  Integrated  Process  Team  (http://www.spawar.naw.mil/).  collected  the 
data  at  10-minute  intervals.  The  data  was  transmitted  to  the  Antarctic 
Meteorological  Research  Center  at  the  University  of  Wisconsin 
(http://amrc.ssec.wisc.edu/)  where  it  was  then  immediately  available  for  our  use. 
The  measurement  of  the  downwelling  solar  radiation  was  by  an  Eppley 
Precision  Spectral  Pyranometer  (PSP)  (The  Eppley  Laboratory,  Inc., 
http://www.eppieviab.ccm/).  which  measures  the  sun  and  sky  irradiance  in  the 
range  of  wavelengths  from  0.285  to  2.8  pm.  The  measurement  of  the 
downwelling  infrared  radiation  was  by  an  Eppley  Precision  Infrared  Radi¬ 
ometer  (PIR),  which  measures  wavelengths  from  3.5  to  50  pm. 

Using  precision  resistors,  we  measured  the  runway  temperatures  at  7 
depths  beneath  the  runway  surface  (4,  6,  8, 10, 12, 14,  and  16  in.).  The 
thermistors  were  placed  at  the  centerline  of  the  runway  and  50  ft  from  the 
edge  of  the  runway  at  3000  and  7000  ft  along  the  runway  for  a  total  of 
four  measurement  locations.  The  thermistors  were  installed  by  first  cut¬ 
ting  a  small  trench  in  the  runway  (Figure  12)  that  extended  from  the  run¬ 
way  centerline  to  the  runway  edge.  Vertical  arrays  of  thermistors  were 
placed  in  the  trench  at  the  runway  centerline  and  50  ft  offset  from  the 
runway  edge.  The  trenches  were  filled  ice  chips  and  snow,  flooded  with 
freshwater,  and  allowed  to  freeze.  Finally,  the  snowcap  was  replaced  on 
the  surface  of  the  ice.  The  thermistors  were  connected  to  a  datalogger  and 
the  information  transmitted  by  satellite  phone  at  1  hr  intervals. 
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Figure  12.  View  of  the  trench  used  for  installing 
temperature  sensors  in  the  runway. 


ERDC/CRREL  TR-15-2 


34 


5  Application 

The  Pegasus  Runway  vertical  temperature  model  was  used  to  simulate  the 
vertical  temperature  profiles  below  the  runway  over  three  consecutive  aus¬ 
tral  summers  during  periods  when  both  meteorological  data  and  runway 
temperature  data  were  available.  Table  3  lists  the  simulation  start  date 
and  time  and  end  date  and  time  of  each  period.  The  model  domain  was 
established  by  dividing  into  layers  the  snow  and  ice  composing  the  runway 
and  the  MIS  beneath  the  runway  down  to  a  15  m  depth.  In  all  simulations, 
the  top  0.5  m  was  divided  into  layers  of  0.01  m  thickness,  the  next  0.5  m 
into  layers  of  0.05  m,  and  the  remaining  depth  into  layers  of  0.10  m.  The 
density,  thermal  conductivity,  and  bulk  extinction  coefficients  were  as¬ 
signed  to  each  layer  based  on  the  vertical  distributions  shown  in  Figures  7, 
8,  and  9.  These  values  were  kept  constant  during  the  simulations  unless 
melting  occurred.  Then  the  density,  thermal  conductivity,  and  heat  capaci¬ 
ty  were  modified  using  equation  (14). 


Table  3.  Simulation  periods. 


Start  Date  Time 

End  Date  Time 

Melt  out 

28  Dec  2011  2200 

29  Jan  2012  1500 

- 

23  Dec  2012  0000 

01  Jan  2013  0000 

29  Dec  2012 

29  Nov  2013  1250 

05  Jan  2014  0320 

31  Dec  2013 

The  first  step  of  each  simulation  was  to  estimate  the  initial  temperature  at 
the  upper  and  lower  boundary  of  each  layer  of  the  solution  grid  at  the  start 
date  and  time  by  using  the  initial  temperature  procedure  described  above. 
The  initial  melt  fraction  was  set  to  zero  for  each  layer  at  the  start  of  each 
simulation.  The  model  was  then  propagated  forward  in  time  using  10- 
minute  time  steps.  The  observed  air  temperature,  wind  speed,  relative 
humidity,  barometric  pressure,  downwelling  shortwave  radiation,  and 
down  welling  long -wave  radiation  were  used  as  inputs  to  the  model.  Dur¬ 
ing  periods  of  missing  data,  the  values  of  these  variables  were  held  con¬ 
stant  using  the  last  available  observation.  The  only  exception  was  the 
downwelling  shortwave  radiation,  which  was  modeled  during  periods  of 
missing  observations  using 


ii  CA 


f  V  \ 


COS  On. 


\r  J 


(30) 
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where 

Ia  =  the  estimated  downwelling  shortwave  radiation; 

r  =  the  daily  distance  between  the  earth  and  sun; 
r0  =  the  mean  distance  between  the  earth  and  sun; 

So  =  the  solar  constant  at  the  mean  Earth-Sun  distance,  r0  = 
1369.3  wm-2; 

6o  =  the  solar  zenith  angle; 

Ci  =  an  empirical  constant  to  account  for  the  impact  of  clouds  (ci  = 
0.75). 

The  value  of  Ci  was  estimated  based  on  the  observations  of  downwelling 
shortwave  radiation  made  at  the  Pegasus  Runway  from  29  October  2010 
until  5  February  2011.  The  solar  zenith  angle  is  a  function  of  the  runway 
location,  the  day  of  the  year,  and  the  time  of  day  and  was  estimated  using 
the  procedure  of  Woolf  (1968). 

Table  4  lists  the  parameters  used  in  the  Pegasus  Runway  temperature 
model.  Only  the  albedo  reduction  factor,  arf,  used  to  reduce  the  surface 
shortwave  albedo  during  times  of  surface  melt  (equation  [6]),  and  ch  the 
empirical  constant  used  to  account  for  the  impact  of  clouds  on  the  esti¬ 
mated  downwelling  solar  radiation  during  periods  of  missing  observa¬ 
tions,  equation  (30),  were  determined  empirically. 


Table  4.  Model  Parameters. 


Parameter 

Value 

Emittance  of  snow 

0.97 

Penetrating  fraction  of  shortwave  radiation 

0.17 

Bulk  extinction  coefficient  blue  ice 

1.5  mr1 

Density  of  ice 

916.8  kg  mr3 

Density  of  water 

1000  kg  mr3 

Thermal  conductivity  of  ice 

2.21  W(m  “C)-1 

Thermal  conductivity  of  water 

0.58  W(m  “C)"1 

Heat  capacity  of  ice 

2114  J  (kg  'C)-1 

Heat  capacity  of  water 

4181 J  (kg  “C)-1 

Latent  heat  of  fusion  of  water 

330000  J  kg-1 

Stefan-Boltzmann  constant 

5.  670  x  IQ-3  W  m-2  k-4 

The  observed  runway  temperatures  were  not  used  in  estimating  the  initial 
vertical  temperature  profile  at  the  start  of  each  simulation  period  or  at  any 
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time  during  a  simulation  period.  The  observed  runway  temperatures  were 
used  only  to  determine  the  performance  of  the  Pegasus  Runway  tempera¬ 
ture  model.  The  results  of  the  simulations  are  described  next. 

5.1  2011-12  season 

The  Pegasus  Runway  temperature  model  was  used  to  simulate  the  2011-12 
season  from  28  December  2011  at  2200  hours  (GMT)  until  29  January 
2012  at  1500  hours  (GMT).  Figure  13  shows  the  observed  meteorological 
data  used  as  input  for  the  model.  There  were  several  periods  of  missing 
data.  The  only  one  of  significant  length  was  from  03  January  2012  at  1540 
hours  until  04  January  2012  at  2000  hours.  During  this  period,  we  esti¬ 
mated  input  parameters  as  described  in  the  introduction  to  this  section. 
Figure  14  displays  the  observed  runway  temperatures,  and  Figure  15  shows 
the  model  results.  In  general,  the  simulation  estimated  the  runway  tem¬ 
peratures  and  their  variations  throughout  the  season  quite  well.  Figure  16 
displays  the  observed  and  simulated  temperatures  at  the  4  in.  depth,  the 
depth  closest  to  the  runway  surface;  and  Figure  17  provides  the  calculated 
temperature  of  the  runway  surface  and  the  calculated  liquid  water  fraction 
of  the  top  surface  layer  and  the  layers  immediately  below  the  top  surface. 
(The  liquid  water  fraction  scale  is  on  the  right  axis  of  Figure  17.  The  liquid 
fraction  result  is  for  the  surface  and  each  layer  immediately  below  the  sur¬ 
face  with  a  liquid  water  fraction  greater  than  zero.)  The  temperature  of 
the  top  surface  reached  o°C  on  several  occasions;  but  the  length  of  time 
that  the  surface  remained  at  o°C  was  relatively  short  for  each  occasion, 
never  lasting  longer  than  a  day.  The  simulation  estimated  that  some  liquid 
water  formed  at  the  surface  and  immediately  below  during  each  event  as 
shown  in  Figure  17. 

It  is  interesting  to  note  that  during  this  season,  the  model  error  goes 
through  three  phases  with  respect  to  time.  These  same  phases  can  be  seen 
at  all  depths,  but  our  description  concentrates  on  the  4  in.  depth  (Figure 
16).  During  the  first  phase,  lasting  from  about  28  December  to  3  January, 
the  model  tended  to  be  too  cold  compared  to  the  observed.  In  the  second 
phase,  from  about  3  January  to  12  January,  the  model  tended  to  be  slightly 
warmer  than  the  observations.  In  the  third  phase,  the  period  after  12  Jan¬ 
uary,  the  model  is  again  slightly  colder  than  the  observations.  The  cause  of 
these  phases  is  not  known.  They  probably  result  from  changes  in  the  run¬ 
way  conditions  that  impacted  the  heat  budget  of  the  runway.  These 
changes  could  have  resulted  from  variations  in  snow  accumulation  on  the 
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runway  through  precipitation,  snow  drift,  or  sublimation;  in  the  compac¬ 
tion,  planning,  or  dragging  operations  performed;  or  in  the  surface  albedo. 


Figure  13.  Meteorological  conditions  during  the  2011-12  season. 
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Figure  14.  Observed  runway  temperatures  during  the  2011-12  season. 


Figure  15.  Modeled  temperatures  in  the  runway  during  the  2011-12  season. 
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Figure  16.  Observed  and  modeled  temperature  at  a  4  in.  depth. 
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Figure  17.  Modeled  temperature  and  liquid  water  fraction  at  the  surface. 
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Table  5  lists  the  error  statistics.  The  simulation  model  results  were  slightly 
negatively  biased  compared  to  the  observed  runway  temperatures. 


Table  5.  Error  statistics  for  the  2011-12  season. 


4  in. 
Depth 

6  in. 
Depth 

8  in. 
Depth 

10  in. 
Depth 

12  in. 
Depth 

14  in. 
Depth 

16  in. 
Depth 

Bias 

-0.85  °C 

-0.67  °C 

-0.54°C 

-0.43  °C 

-0.35  °C 

-0.17  °C 

-0.02  °C 

StDev 

0.69  °C 

0.61°C 

0.54  °C 

0.48  °C 

0.43  °C 

0.39  °C 

0.37  °C 

5.2  2012-13  season 

The  Pegasus  Runway  temperature  model  was  used  to  simulate  the  2012- 
13  season  from  23  December  2012  at  0000  hours  (GMT)  until  05  January 
2014  at  0320  hours  (GMT),  when  the  runway  melted  to  the  point  that  the 
top-most  runway  temperature  sensors  were  exposed  to  air  and,  therefore, 
the  data  was  no  longer  valid.  Figure  18  shows  the  observed  meteorological 
data  used  as  input  for  the  model.  There  were  no  significant  periods  of 
missing  data.  Figure  19  displays  the  observed  runway  temperatures.  Note 
that  there  is  a  general  decline  in  the  observed  runway  temperatures  start¬ 
ing  from  their  first  reading  until  approximately  23  December.  This  two- 
day  decline  is  thought  to  be  an  artifact  resulting  from  the  procedure  used 
to  install  the  temperature  probes  into  the  runway.  It  can  also  be  seen  that 
the  observed  runway  temperatures  increase  more  or  less  continuously  un¬ 
til  01  January  when  the  Pegasus  Runway  was  closed  to  operations. 

Figure  20  displays  the  model  results.  In  general,  the  simulation  estimates 
of  the  runway  temperatures  are  reasonable.  The  model  was  able  to  repro¬ 
duce  quite  well  the  temperature  rise  in  the  runway,  and  its  timing,  that  led 
to  the  runway  shutdown.  Figure  21  displays  the  observed  and  simulated 
temperatures  at  the  4  in.  depth,  the  depth  closest  to  the  runway  surface. 
Figure  22  displays  the  calculated  temperature  of  the  runway  surface  and 
the  calculated  liquid  water  fraction  of  the  top  surface  layer  and  the  layers 
immediately  below  the  top  surface.  (The  liquid  water  fraction  scale  is  on 
the  right  axis  of  Figure  22.  The  liquid  water  fraction  result  is  for  the  sur¬ 
face  and  each  layer  immediately  below  the  surface  with  a  liquid  water  frac¬ 
tion  greater  than  zero.)  The  temperature  of  the  top  surface  reached  o°C  on 
27  December  and  consistently  reached  o°C  on  each  of  the  following  days. 
Some  liquid  water  was  estimated  to  form  at  the  surface  and  immediately 
below  during  each  event  as  shown  in  Figure  22.  Table  6  lists  the  error  sta- 
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tistics.  The  simulation  results  were  negatively  biased  compared  to  the  ob¬ 
served  runway  temperatures. 


Table  6.  Error  statistics  for  the  2012-13  season. 


4  in. 
Depth 

6  in. 
Depth 

8  in. 
Depth 

10  in. 
Depth 

12  in. 
Depth 

14  in. 
Depth 

16  in. 
Depth 

Bias 

-1.09  °C 

-1.01°C 

-0.83  °C 

-0.70°C 

-0.53  °C 

-0.38°C 

-0.27  °C 

StDev 

1.29  °C 

0.99  °C 

0.76°C 

0.62  °C 

0.58°C 

0.59  °C 

0.63  °C 
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Figure  19.  Observed  temperatures  in  the  runway  during  the  2012-13  season. 
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Figure  21.  Observed  and  modeled  temperature  at  a  4  in.  depth. 


Figure  22.  Modeled  temperature  and  liquid  water  fraction  at  the  surface. 
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5.3  2013-14  season 

The  Pegasus  Runway  temperature  model  simulated  the  2013-14  season 
from  29  November  2013  at  1250  hours  (GMT)  until  1  January  2013  at 
0000  hours  (GMT),  when  the  runway  melted  to  the  point  that  the  topmost 
runway  temperature  sensors  were  exposed  to  air  and,  therefore,  the  data 
was  no  longer  valid.  Figure  23  shows  the  observed  meteorological  data 
used  as  input  for  the  model.  There  were  a  number  of  periods  of  missing 
data,  including  one  significant  four-day  period  from  20  January  until  23 
January.  In  addition,  we  determined  that  the  observations  of  the 
downwelling  solar  radiation  were  in  error;  therefore,  we  used  estimated 
values  for  the  entire  simulation. 


Figure  23.  Meteorological  conditions  during  the  2013-14  season. 
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Figure  24  shows  the  observed  runway  temperatures,  which  begin  on  6  De¬ 
cember  2013.  As  in  the  previous  season,  there  is  a  general  decline  in  the 
observed  runway  temperatures,  starting  from  their  first  reading  and  last¬ 
ing  approximately  one  day.  As  before,  we  assume  this  to  be  an  artifact  re¬ 
sulting  from  the  procedure  used  to  install  the  temperature  probes  into  the 
runway.  Additionally,  the  observed  runway  temperatures  increase  starting 
at  the  end  of  December  and  continuing  until  04  January  when  the  Pegasus 
Runway  was  closed  to  operations.  Figure  25  displays  the  model  results. 
The  model  simulation  period  started  on  29  November,  prior  to  the  start  of 
the  observations.  In  general,  the  simulation’s  estimate  of  the  runway  tem¬ 
peratures  is  very  good,  especially  considering  the  significant  periods  of 
missing  meteorological  data.  The  model  was  able  to  effectively  reproduce 
the  runway’s  temperature  rise  (and  its  timing)  that  lead  to  its  shutdown. 
Figure  26  shows  the  observed  and  simulated  temperatures  at  the  4  in. 
depth,  the  depth  closest  to  the  runway  surface.  Figure  26  provides  the  cal¬ 
culated  temperature  of  the  runway  surface  and  the  calculated  liquid  water 
fraction  of  the  top  surface  layer  and  the  layers  immediately  below  the  top 
surface.  (The  liquid  water  fraction  scale  is  on  the  right  axis  of  Figure  26. 
The  liquid  water  fraction  result  is  for  the  surface  and  each  layer  immedi¬ 
ately  below  the  surface  with  a  liquid  water  fraction  greater  than  zero.)  The 
temperature  of  the  top  surface  reached  o°C  on  several  occasions  and  then 
was  consistently  at  o°C  from  30  December  onwards.  Liquid  water  was  es¬ 
timated  to  form  at  the  surface  and  immediately  below  during  each  occa¬ 
sion  as  shown  in  Figure  27.  Table  7  lists  the  error  statistics.  The  simula¬ 
tion  results  were  slightly  negatively  biased  above  the  8  in.  depth  compared 
to  the  observed  runway  temperatures  and  slightly  positively  biased  below 
that  depth. 


Table  7.  Error  statistics  for  the  2013-14  season. 


4  in. 
Depth 

6  in. 
Depth 

8  in. 
Depth 

10  in. 
Depth 

12  in. 
Depth 

14  in. 
Depth 

16  in. 
Depth 

Bias 

-0.53  °C 

-0.19  °C 

-0.08°C 

0.13  °C 

0.32  °C 

0.49  °C 

0.72  °C 

StDev 

0.79  °C 

0.71°C 

0.62°C 

0.55°C 

0.51°C 

0.47  °C 

0.44  °C 
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Figure  24.  Observed  temperatures  in  runway  during  2013-14  season. 


Figure  25.  Modeled  temperatures  in  runway  during  2013-14  season. 
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Figure  26.  Observed  and  modeled  temperature  at  a  4  in.  depth. 
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Figure  27.  Modeled  temperature  and  liquid  water  fraction  at  the  surface. 
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6  Summary 

The  Pegasus  Runway,  located  on  the  McMurdo  Ice  Shelf  (MIS)  13  km 
south  of  McMurdo  Station,  is  a  unique  runway  composed  of  a  compacted 
snowcap  supported  by  a  freshwater  ice  layer  approximately  4.7  m  thick. 
The  runway  is  susceptible  to  weakening  and  reduction  in  the  snowcap 
strength  caused  by  warm  weather  and  solar  radiation.  In  the  present 
work,  we  describe  the  development  of  a  temperature  model  of  the  upper 
layers  of  the  Pegasus  Runway.  This  model  is  a  major  component  of  the 
system  to  forecast  the  occurrence  of  strength  reduction  in  the  airfield, 
which  will  be  described  in  following  reports. 

To  model  the  temperatures  in  the  Pegasus  runway,  it  is  necessary  to  have  a 
description  of  density,  thermal  conductivity,  specific  heat,  and  the  broad¬ 
band  shortwave  radiation  extinction  of  the  snow  and  ice  that  compose  and 
underlay  the  runway.  The  surface  albedo  must  be  known,  also.  Measure¬ 
ments  of  these  physical  properties  in  the  field  are  limited  to  density  along 
with  a  few  observations  of  broadband  shortwave  radiation  extinction.  The 
direct  measurements  of  density  made  at  the  runway  by  us  and  others  were 
confined  to  the  top  5  m.  Based  on  measurements  made  by  others  in  the 
vicinity  of  the  Pegasus  Runway,  we  inferred  the  density  deeper  than  5  m. 
We  then  estimated  the  other  physical  properties  based  on  the  measured 
and  estimated  density  and  the  layer  type. 

In  all,  there  are  four  distinct  layers  of  snow  and  ice  at  the  runway  between 
the  surface  and  the  bottom  of  the  MIS:  snowcap,  superimposed  ice,  firn, 
and  brine-saturated  firn.  Of  the  four  distinct  layers  of  snow  and  ice,  the 
brine-saturated  firn  presented  the  greatest  challenge  with  regard  to  heat 
transfer  calculations.  The  presence  of  brine  has  a  dramatic  impact  on  the 
volumetric  specific  heat,  especially  at  temperatures  close  to  the  ice/sea  wa¬ 
ter  equilibrium  temperature  and  high  salinities.  In  practice,  the  specific 
heat  of  brine-saturated  firn  could  be  estimated  if  the  brine  fraction  and 
brine  salinity  were  known.  These  properties  are  not  known  beneath  the 
Pegasus  Runway.  If  we  assumed  that  the  brine-saturated  firn  was  formed 
when  the  original  dry  firn  (with  an  estimated  density  determined  either  by 
depth  or  overburden  pressure)  was  flooded  with  seawater  (with  known  sa¬ 
linity),  we  could  estimate  the  properties  of  the  brine-saturated  firn.  Un¬ 
fortunately,  measurements  of  the  properties  of  the  brine-saturated  firn  in 
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the  vicinity  of  Pegasus  Runway  do  not  present  a  consistent  picture  of  the 
original  brine  salinity.  We  approximate  the  thermal  properties  of  the 
brine-saturated  firn  based  on  the  properties  of  ice  with  bubbles  as  de¬ 
scribed  by  equation  (8).  This  approximation  is  justified  given  the  relative 
depth  of  the  brine-saturated  firn.  Below  15  m  to  the  base  of  the  MIS,  the 
temperature  does  not  change  with  time;  and  the  specific  heat  does  not  ap¬ 
pear  in  the  thermal  calculations.  Between  15  m  and  the  top  of  the  brine- 
saturated  firn  at  9  m,  the  ice  is  cold;  and  the  brine  has  less  impact  on  the 
thermal  properties.  Overall  we  assumed  that  the  saturated  brine  layer  has 
the  thermal  conductivity  of  ice  with  a  density  of  880  kg  m~3  and  a  result¬ 
ant  thermal  conductivity  of  2.13  W  (m  °C)_1. 

The  estimated  thermal  properties  of  the  four  distinct  layers  of  snow  and 
ice  at  the  runway  allowed  the  model  to  accurately  simulate  the  observed 
runway  temperatures,  showing  that  the  assumptions  made  for  the  proper¬ 
ties  of  the  brine-saturated  firn  were  appropriate.  Our  technique  for  esti¬ 
mating  the  initial  temperature  profile  at  startup  worked  well  by  using  em¬ 
pirical  modifications  to  the  classic  mathematical  description  of 
propagation  of  periodic  oscillations  of  surface  temperature.  The  broad¬ 
band  shortwave  albedo  of  the  runway  was  developed  based  on  field  meas¬ 
urements  made  during  the  winter  of  2010-11  of  the  downwelling  and 
upwelling  shortwave  radiation.  We  found  that  the  parameterization  pro¬ 
posed  by  Dickinson  et  al.  (1986)  described  the  albedo  well  as  a  function  of 
the  solar  angle.  We  reduced  the  surface  albedo  at  a  time-constant  rate 
when  the  liquid  water  fraction  of  the  surface  layer  was  greater  than  zero. 
This  allowed  accurate  simulation  of  the  meltout  of  the  runway  in  the 
2012-13  and  2013-14  seasons. 

The  model  of  the  temperature  within  the  Pegasus  Runway  compares  very 
well  to  the  temperatures  measured  in  the  field.  In  almost  all  cases,  the  av¬ 
erage  error  at  each  depth  is  less  than  i.o°C.  In  general,  the  average  error 
(bias)  becomes  less  negative  with  depth.  Of  the  three  simulated  seasons, 
the  2011-12  season  provides  the  longest  period  of  time  for  comparison.  It 
is  interesting  to  note  that  during  this  season,  the  model  error  goes  through 
three  phases  with  respect  to  time.  These  same  phases  can  be  seen  at  all 
depths.  During  the  first  phase,  lasting  from  about  28  December  to  3  Jan¬ 
uary,  the  model  tended  to  be  too  cold  compared  to  the  observed.  In  the 
second  phase,  from  about  3  January  to  12  January,  the  model  tended  to 
slightly  warmer  than  the  observations.  In  the  third  phase,  the  period  after 
12  January,  the  model  was  again  slightly  colder  than  the  observations.  We 


ERDC/CRREL  TR-15-2 


50 


do  not  know  the  cause  of  these  phases,  but  they  probably  result  from 
changes  in  the  runway  conditions  that  impact  the  heat  budget  of  the  run¬ 
way  but  that  are  not  adequately  captured  in  the  model.  These  changes 
could  have  resulted  from  changes  in  snow  accumulation  on  the  runway 
through  precipitation,  snow  drift,  or  sublimation;  changes  in  the  compac¬ 
tion,  planning,  or  dragging  operations  performed;  or  changes  in  the  sur¬ 
face  albedo. 
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